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■ We compute the inclusive differential cross section production of the pseudo-scalar 

meson rf in high-energy proton-proton (pp) and proton-nucleus (pA) collisions. We 
use an effective coupling between gluons and rj meson to derive a reduction formula 



that relates the rj' production to a field-strength tensor correlator. For pA collisions 
we take into account saturation effects on the nucleus side by using the Color Glass 
Condensate formalism to evaluate this correlator. We derive new results for Wil- 
Q h i son line - color charges correlators in the McLerran-Venugopalan model needed in 

Q-ij the computation of rj production. The unintegrated parton distribution functions 

^ . are used to characterize the gluon distribution inside protons. We show that in pp 

collisions, the cross section depends on the parametrization of unintegrated parton 
distribution functions and thus, it can be used to put constraints on these distri- 

-xj- : 

, butions. We also demonstrate that in pA collisions, the cross section is sensitive to 

fN) ' saturation effects so it can be utilized to estimate the value of the saturation scale. 
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I. INTRODUCTION 

• rH 

X 

One of the main challenges in high energy hadronic collisions is the understanding of 
particle production. The theoretical description of these complex phenomena involves both 
many-body physics and the theory of strong interactions. There are many approaches that 
aim toward a better comprehension of these topics. One of the most successful technique 
is perturbative quantum chromodynamics (pQCD) where one studies the limit where the 
coupling constant is small and where the usual loop expansion can be used in principle. 
Because QCD is an asymptotically free theory, this happens when the exchanged momenta 
are large compared to the QCD scale Aqcd- Even in that regime however, physical observ- 
ables computed using this machinery suffer from infrared divergences that spoil the naive 
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perturbative expansion. These have to be resummed and this leads to various factorization 
formalisms like the collinear factorization and the /c^-factorization. In collinear factoriza- 
tion, meaningful physical quantities are obtained in terms of parton distribution functions 
(PDF). These distributions characterize the non-perturbative (large distance) physics and 
have to be determined experimentally from a fit to structure functions. This formalism 
can be applied on system where the typical exchanged momentum Q is hard, which means 
that it satisfies the inequality given by Aq CD <C Q 2 ~ s where Aq C d ~ 0.2 GeV is the 
QCD scale and y/s is the center of mass energy. Note that Q 2 ~ M\ = M 2 + k\, which 
means that unless the mass of the produced particle M is of the order of yfs, the collinear 
factorization is valid only for very large kj_. This can be relaxed in the fc^-factorization for- 
malism which considers semihard collisions, meaning that the typical exchanged momentum 
obeys Aq CD C Q 2 C s. The resummation implemented in this approach takes care of large 
contributions that look like ^]ji(Q 2 / AQ CD )a s ] , [\n(Q 2 /Aq CD ) ln(l/a;)a s ] " and [ln(l/x)a; s ] ra 
This technique is used successfully to compute the production of many kinds 
of particles in high energy proton-proton collisions like heavy quarks and in a 

number of other processes (see 0, |l0, U| for reviews of many applications). It is also used 



for predictions of Higgs boson production 12j, |l3j]. In this article, we use this formalism to 



compute the rf inclusive cross section at the RHIC energy. 

When there is a nucleus involved in a collision at very high energy, there are new effects 
not included in the previous approaches due to the high density of gluons resulting from 
the emission enhancement at small- a; (where x is the momentum fraction). These effects 
introduce a new scale Q s called the saturation scale at which the probability of having 
interactions between gluons of the same nucleus becomes important. At this transverse 
momentum scale, the gluons recombine and this slows down the growth of partons at smaller 
x. A naive estimation of Q s shows that it depends on x and the number of nucleons A like 
Q 2 ~ A s x~ x Jj], [l5] so at small enough x or large enough A, the saturation scale is hard 
(Q 2 ^> Aq CD ). When the typical exchanged momentum is smaller than the saturation 
scale such as Aq CD <C Q 2 < Q 2 <C s, saturation effects have to be taken into account 
even if the system is still in the perturbative regime. This can be achieved in the Color 
Glass Condensate (CGC) formalism which is a semi-classical effective theory where the non- 
linearities are dealt with by solving exactly the Yang-Mills equation of motion. This takes 
care of gluon recombinations and introduces the effects of saturation in observables. 

In this article, we are using the CGC to compute the inclusive differential cross section of 
rj meson in pA collisions at the RHIC energy (y/s = 200 GeV). The main goal of this work 
is to look at the effect of saturation in 77' production to validate the CGC approach and 
estimate the value of the saturation scale by comparing our predictions with experimental 
data. The rj is a pseudoscalar meson with a mass of M = 0.957 GeV, a decay width of 
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0.203 MeV and quantum numbers of I G (J PC ) = + (0 
features of rf is that it couples to the QCD anomaly 17 



|lq| . One of the most important 
18]. One way to implement and 
model this physics is by introducing an effective interaction between gluons and rf mesons. 
This was done in 19] , where the authors are proposing a vertex that couples two gluons and a 
rf meson (ggrf) to explain B mesons decay (B — > rf +X). This vertex contains a form factor 
that depends generally on gluons and rf momenta and that can be related to the rf wave- 
function. The structure of this vertex was investigated thoroughly using various techniques 



like the hard scattering and the running coupling approaches. |20L |21| . [22j, [23|, |24], |25( . We 
use these results on the gluons-?/ coupling to study the rf production mechanism based on 
gluon fusion. 

The rf production in pp collisions at high energy was studied in 26] where the exclusive 
cross section for the diffractive process p + p — > p + p + rf is computed. In our study, we 
focus on the inclusive production mechanism p + p — ■> rf + X which shares similar features 
with this previous analysis. The first attempt to compute rf production in high energy 
pA collisions was done by one of the present author in {27]]. In this study, the collinear 
factorization is used to compute the cross section at RHIC by including intrinsic transverse 
momentum in the PDF with a Gaussian distribution. Based on physical arguments, the 
width of the Gaussian, which represents the typical transverse momentum of gluons inside 
the nucleus, is chosen to be Q s . The authors show that the rf production is sensitive to the 
saturation scale implemented in this way. However, they acknowledge that their calculation 
can be improved because they use the collinear formalism outside its range of validity. The 
goal of this article is to revisit the rf production with a more rigorous approach by doing a 
full CGC computation that includes recombination effects more realistically. 

The computation of meson production in the CGC was undertaken in the past using 
mostly an "hybrid" approach where the proton and the nucleus are described by the collinear 
factorization and the CGC respectively. In this formalism based on pQCD-like techniques, 
the fragmentation function of collinear factorization is convoluted with the gluon or quark 
cross section computed in the CGC formalism. This is suitable for well-known mesons like 

gmenta- 



ra. 
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30] 



pions for which a wealth of experimental data have been measured and for which 
tion functions are well-known. Pion production for pA collisions is computed in 
using this methodology. Contrary to pions, the data in high energy hadronic collisions for rf 
is scarce, so another approach is required to take care of hadronization effects and internal 
structure of the rf meson. In 31( , an effective theory is used to estimate the tensor meson 
production in pp collisions. We use a similar approach in this article where the interaction 
between gluons and rf is described by an effective theory. As discussed previously, we in- 
clude these effects in an effective vertex that includes a form factor. As will be shown in 
this article, this can be implemented easily in the CGC formalism. 
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We consider only the case of pp and pA collisions at RHIC. For nucleus-nucleus (^4^4) 
collisions, the total number of //-mesons produced by semihard collisions in the first instants 
(for t < 1 fm/c) should be important and the saturation effects would also be present. 
Experimentally however, they cannot be detected because most of them decay inside the 
medium created by the collision. This is because the rf mean lifetime, which is about 
t v r pa 4.93 fm/c, is smaller than the time where the medium exists, which is from 1 fm/c up 
to 10 fm/c. Moreover, by considering pp and pA, we avoid all the complications that would 
result from the creation of the medium which include the understanding and modelling of 
the quark-gluon plasma properties. Finally, there are analytical solutions for the gauge field 
in pp and pA collisions, while the analytical solution in AA is still elusive. For these reasons, 
our present analysis is only applied to pp and pA collisions. 

This article is organized as follows. In section [TT] we describe the effective vertex used 
throughout the rest of the article. In section IHH we show how to compute rf production in 
the CGC for pA and pp collisions. We start by deriving a reduction formula that relates the 
cross section to a correlator of field-strength tensors. This correlator is then evaluated to 
leading order in pA collisions and the result can be interpreted in terms of physical processes. 
We also show how the fcj_-factorized cross section for pp collisions can be recovered in the low 
density limit of pA cross section. In section IIV} we compute the correlators appearing in the 
expression of the pA cross section using the McLerran-Venugopalan model. In section |V] we 
evaluate numerically the cross section for pA and pp and discuss the range of validity of our 
computation. Sections III and IV contain a lot of technical details. The reader interested 
in results can jump directly to section V. 

Throughout the article, we use both light-cone coordinates defined by 

, p° + p 3 _ p° — p 3 
v = — t~ ; v = — T~ I 1 ) 

y/2 V2 

and Minkowski coordinates. It should be clear by the context which one is used. We also 
use the metric convention = (1,-1,— 1,-1). 



II. EFFECTIVE THEORY 



The effective theory used in this article couples gluons and the rf meson. In momentum 
space, the g*g*rj' effective vertex (where g* means off-shell gluon) is given by 

V^(M,p, q) = iF(p 2 , q 2 , M 2 )5 ab e^ a p p q a (2) 

where e^ upa is the Levi-Civitta antisymmetric tensor, M is the rf mass, a and b are color 
indices, p p and q a are gluon momenta and F(p 2 , q 2 , M 2 ) is the rf form factor. The explicit 
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expression of the interaction vertex have been studied in a number of articles where different 
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251 ] . To get a 



parametrizations of the form factor can be found 
first approximation of rf production and because we are mostly interested in making a 
comparative study between pp and pA collisions, we use a simple expression given by 26 ] 

M 4 



F(p 2 ,q z ,M 2 ) = H 



(3) 



(M 2 — p 2 )(M 2 — q 2 ) 

where Hq = F(0, 0, M 2 ). To get a better approximations of rf production, other parametriza- 
tions should be used. In the limit of no gluon virtualities (p 2 , q 2 



0), the form factor 
It is given by 

H Q = F(0, 0, M 2 

The main applications of this coupling are related mostly to B and T decay where 



is a constant that can be fixed by looking at the decay of ip — > rf + 7 
r2 ^ ~ 1.8 GeV 0. 



processes such as g* — > g + rf and g* + g* — > rf are considered [l9|, |32j, |33, 



34 



35j. More 



rece ntly , gluon fusion was used to compute rf production in high energy hadronic collisions 
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361 ] and from a thermalized medium 371 ] . 



It is convenient for our purpose to consider the interaction Lagrangian given by 
C iat (x) =\f d A yd A zF [{x - y)\ (x - z) 2 ,M 2 } G^(y)G a ^(z) v (x). 



(4) 



that reproduces the vertex Eq. (j2J) in the perturbative expansion. As seen in the next 
section, this can then be used to derive a reduction formula. Here, G£ u (x) is the usual 
field-strength tensor given by 



0?{x) = d»Al{x) - (TA»{x) - gf abc A»{x)A v c {x) 



(5) 



where A% is the gauge field of gluons and G^ u (x) = ^ vpa G a ^ a [x) is the dual field-strength 
tensor. The Lagrangian is non-local because the vertex includes a form factor. It can be 
easily seen that C mt {x) contains three types of vertices, namely ggr]', gggrf and ggggrf . At 
leading order however, only the first one is necessary and considered in this article. 



III. PRODUCTION OF rf FROM THE CGC 

In collisions at very high energy, the wave function of nuclei is dominated by gluons 
that have small longitudinal momenta (soft gluons) because of the emission enhancement 
at small-x. The CGC is a semi-classical formalism that describes the dynamics of these 
degrees of freedom. In this approach, the hard partons, which carry most of the longitudinal 
momentum, and soft gluons which have small longitudinal components, are treated differ- 
ently. Because the occupation number of the soft gluons is large, classical field equations 
can be employed to understand their dynamics. The hard partons act as sources for these 
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d and are no longer interacting with the rest of the system (for reviews of CGC, 



aal) 



In this formalism, computing a physical quantity involves two main steps. The first one 
is to solve the Yang-Mills equation of motion 



[D ll ,F^{x)] = J v [ 



x 



(6) 



where the current J%(x) = 5 u+ 5(x )p p n {x \ ) + 5 U S(x + )pA } a( x ±) represents random static 



sources localized on the light-cone 14], Ilia] and _D M = — igA^ is the covariant derivative. 



The functions p Pi a(x±) are color charge densities in the transverse plane of the proton and 
nucleus respectively. The next step is to take the average over the distribution of color 
charge densities in the nuclei with weight functionals W Pi a[Pp,a]- For any operator that can 
be related to color charge densities, this can be written as 



(O) 



Vp p Vp A 0[pp, pa]W p [x p , p p ]W a [xa, Pa]- 



(7) 



Computing the weight functional is a highly non-perturbative procedure so it usually involves 
approximations based on physical modelling. In the limit of a large nuclei at not too small x, 
it can be approximated by the McLerran-Venugqpalan (MV) model, which assumes that the 



partons are independent sources of color charge [39|, |40j. Within this assumption, the weight 
functional W ^I/jJ is a xa independent Gaussian distribution and the two point correlator 
is simply QHS&I 



(pA, a ( x ±)PA,b(y±)} = 5 ab p A d ( x ± - y± 



(8) 



where p 2 A = A/2ttR 2 is the average color charge density and R is the radius of the nucleus. 
It is assumed here that the nucleus has an infinite transverse extent with a constant charge 
distribution. Edge effects can be included by changing p 2 — > p 2 (x±) and by choosing a suit- 
able transverse profile. Throughout this article, we only consider the constant distribution 
case. 

Within the MV model, the weight functional does not depend on longitudinal coordinates 
and therefore, the model is boost invariant. This however can be relaxed by considering 
the quantum version of the CGC. In that theory, quantum radiative corrections become 
important below a certain scale xq ~ 0.01. These corrections can be resummed b y using 
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42 



a renormalization group technique which leads to the JIMWLK evolution equation 
43l . |44], |45(. In the quantum CGC, the weight functionals Wi^ [pi.,2] obey this non-linear 
evolution equation in x. Because the MV model is valid in the range x ~ 0.01 — 0.1, it 
can be used as an initial condition for the evolution at smaller x. In this article however, 
we consider only the regime where the MV model is valid and do not consider the small- x 
evolution although it could be done in principle. 
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On the proton side, the average computed with W p [x p , p p ] can be related to the uninte- 
grated parton distribution function (uPDF) <\>\ like 

P± + q± 



x / ^ ± ^^x j^gHgj (9) 



2 d0i 

y± ™ = <w,p±) (io) 



where iV c is the number of color. By construction, the uPDF obeys 

d 2 y. 

and is normalized such that 

.2 



01,2(^,P±) ~ ^(X,^ 2 ) (11) 







where xG(x, /i 2 ) is the collinear parton distribution function and /x 2 is the factorization scale. 
The uPDF can be obtained from a fit to structure function and evolved to the desired value 
of x p , Q 2 and p 2 ± using evolution equations such as the BFKL or the CCFM equations. 

One important ingredient is missing for the computation of rf meson production cross 
section. We need a relation between the cross section and a correlator that can be evaluated 
using the CGC formalism. This is done in the next section using a reduction formula and 
the effective theory. 

A. Reduction Formula and the Cross Section 

The computation of rf mesons from the CGC can be calculated from a reduction formula. 
The starting point is the expression of the average number of rf produced per collisions given 
by n = Y^nLi n Pn where P n is the probability to produce n particles. This can be converted 
to an equation in terms of creation/annihilation operators that can be evaluated in quantum 



field theory. This is given by 46] 



(2vr) 3 2£ fc J| = (0 in | at ut (k)a out (k) |0 in ) (12) 
where |0 in ) is the in vacuum. Then, the standard LSZ procedure can be used to write this 



as 47] 



(2-nf2E k — = -J d 4 xd A ye lk ^ [d 2 x + M 2 ] [d 2 + M 2 ] (0 in | fj(x)rj(y) |0 in ) (13) 

where Z is the wave function normalization and where we assumed the asymptotic conditions 
limt^i-oo fj(x) = VZfj out ; m (x) for the rf field operator in Heisenberg representation. It is 
possible to use the equation of motion of r){x) given simply by 

(0 2 + M 2 )fj{x) = f{x) (14) 
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where we defined f(x) = ~ J d A yd A zF [(x — y) 2 , (x — z) 2 , M 2 } G% u (y)G a>fll/ (z) to rewrite the 
reduction formula in a convenient way. We get finally that 

(i 71 a a 

where T(/c) is the Fourier transform of T(x) evaluated at a rj' meson on-shell momentum 
and where we set Z = 1. The angular brackets (O) here indicates expectation value of O in 
the initial state. 

The only assumptions used in deriving Eq. (|15p are: 

• There are no if mesons in the initial state. 

• The rf is produced on-shell. 

The first assumption is justified by the fact that in high-energy collisions, the number of 
rf in a hadron before the collision (in the initial state) is negligible. This allows us to 
use the in vacuum and the fact that a\ n (k)\Q\ n ) = to simplify the reduction formula. 
Using the second assumption, we can treat the r( meson as a stable particle which can 
be produced on-shell and which is well described by the free spectral density that looks 
like p(M 2 ) ~ 5(p 2 — M 2 ). Therefore, by making this assumption, it is possible to use the 
asymptotic conditions described earlier. However, 7/-mesons are resonances, so the spectral 
density should look rather as a Breit-Wigner function p(M 2 ) ~ T/[(p 2 — M 2 ) 2 + M 2 Y 2 ] 
where V is the decay width. These effects however are taken into account by the form factor 
F(q,p,M). 

Eq. (II 5p is the main result of this section. It relates the average number of mesons 
produced to a correlator of field strength tensors. This correlator can then be evaluated using 
any analytical or numerical methods. The average (...) depends on the system studied. 
Looking at a plasma at equilibrium, it could be computed using finite temperature field 
theory or the AdS/CFT correspondence. These two formalisms are relevant to nucleus- 
nucleus collisions where a medium at equilibrium is created. We are interested here in pA 
and pp collisions where no such medium is formed so these techniques are not pursued in 
this study. Rather, we use the CGC which describes initial state and saturation effects in 
high-energy hadronic collisions. 

Having expressed the average number of rj' produced in terms of a correlator of field 
strength tensors, it is possible to compute the inclusive cross section in the CGC formalism 
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which is given by [46|, |48 ] 



(2w) 3 2E k 



da 
~¥k 



d 2 b ± (2n) 3 2E, 



dn(b± 
: d 3 k 



d 2 b± J Vp p Vp A T*[p p , p A ]T[p p , p A ] 
xW p [x p , p p ]WA[pA;bx\- 



(16) 



In this expression, b± is the impact parameter. The fields T are functionals of the source 
once the Yang-Mills equation of motion of the gauge field is solved (see Eq. (JHJ) for the 
expression of the field strength tensor as a function of the gauge field). 



B. Cross section in pA Collisions 



In pA collisions, there are two saturation scales (one for the proton (Q p ) and one for 
the nucleus (Qa)) that satisfy Q p < Qa- When the transverse momentum of the rj' is 
small enough, the nucleus is in a saturation state while the proton is not because we have 
Qp < Q 2 = M\ < Q 2 A (remember that M\ = M 2 + k\ is the transverse mass of the 
T]'). In that case, the system is semi-dilute, meaning that one of the source is strong (or 
equivalently, the typical transverse momentum is small) and obeys PA,a(k±)/k 2 _ ~ 1 while 
the other source is still weak p P! a(k±)/kj_ <C 1 0, 4^- The weak source can be used as 
a small parameter to solve the Yang-Mills equation perturbatively. The solution of the 

y to all orders in PA.a(k±)/k 2 _ and to first order in 



52 



53]. We use here the solution in the light-cone 



gauge field can be computed analytical 
Pp,a(k±)/k 2 _ in different gauges [50l. l5ll. 
gauge of the proton 52J but in Appendix [Cj we perform the same calculation in covariant 
gauge to show that our result is gauge invariant. 



1. Gauge Field and Power Counting 

In the light-cone gauge of the proton (A + = 0) with a nucleus in covariant gauge moving 
in the negative z direction, the solution of the gauge field in pA collisions can be separated 
in three parts A^{k) = A^^k) + A^ a (k) + A^ A where A£ a (k) is the field associated with 
the proton (of 0(p p )), A Aa (k) is the field associated with the nuclei (of O(pa) ~ 0(1)) and 
Ap Aa (k) is the field produced by the collision (of 0(p p p A ) ~ 0(p p )) [521 ] . Note that the field 
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A A is strong and satisfies Ap 1 , A^ A <C A A . The explicit solution is given by [52l. \54 \ 



A p>a (k) = togS{k-) J* ^ (17) 



A~ Aa (k) = 2ng5(k+)^j^ (18) 
ig f d 2 p ± 



A p a ^ k 2 + ik+e J (2tt) 2 



k % ^ V 1 

(k+ + ie)(k~ + ie) p\ 



Pp,&(Pj 



x [U ab (k ± -p ± ) - (27i) 2 5 2 (k ± - P± )5 ab ] (19) 

where U ab (k±) is a Wilson line in adjoint representation defined in Eq. fl35|) . g is the usual 
QCD coupling constant, 5 ab is the Kronecker delta in color space and f abc is the structure 
constant of the SU(N C ) group. The component A~ A a {k) is non-zero and is related to A l pA a (k) 
but it does not appear in the final expression of the cross section so it is not needed in the 
computation of rf production. All the other components are zero. 

The production cross section of rf mesons is related to a field strength tensor correlator 
given by 

B(k) = J d ^F{p 2 y 2 )F\ q 2 ,ql)(Gl u MGT{k (20) 

where P2,<i2 = k — p,q. In principle, a correlator like this contains contributions from all 
orders in both sources. Because the proton source is weak, it is possible to simplify this 
considerably using a power counting argument to isolate the leading order contribution. 
For the sake of this power counting argument, we use a? s which denotes terms having n 
gauge fields A^ and where ~ O(pa) and a p , cl p a ~ 0(p p ). At first, let us consider only 
the contributions from the abelian part of the field-strength tensor. The terms in these 
contributions have four powers of gauge field such as -B a beiian ~ a tot ~ ( a A + a P + (i p a) 4 - 
Naively, one would expect the dominant contribution to come from terms that have many 
powers of the nucleus gauge field like a A ~ 0(p A ) and a 3 A a p ,a A a p A ~ 0(p A p p ). However, 
these terms vanish because of the Lorentz structure. For example, a typical term in the 
abelian contribution would look like T abelian ~ e^ upa p fi q u A p A' J . When we sum on indices, 
this kind of term will contain at most one strong gauge field A~ ~ a a- Thus, the dominant 
contributions are like a 2 A a 2 ,a A a 2 A ~ 0{p\p 2 ). Using a similar argument, it is possible to 
show that the non-abelian part have no leading order contribution in the sense that it is 
at least -B n on-abeiian ~ 0(.PaPp) ^ O(PaPp)- The possible higher order contributions like 
a A a p ~ O(PaPp) f° r example also vanish because of the Lorentz structure of the correlator. 
This is because the typical non-abelian contributions look like 7^ on _ abelian ~ e tlupa p tJ, A u A p A a 
and ^non-abciian ~ ^ i ivpaA p A v A p A a . Once the Lorentz indices are summed, the second typical 
term T" = because in this gauge, A + = 0. For T', it contains only one strong field A~ ~ a p A 
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but it contains two powers of weak field like A % ~ a p , cl p a- Thus, when it is squared, it gives 
at most a contribution of J B non _ abclian ~ 0{p\p z p ). 

C. Evaluation of the Correlator 

Using the explicit expression of the field strength tensor in terms of gauge field and 
keeping only the dominant and non-zero contributions, the correlator can be written as 

B(k) = 64 1 ^^F(-pi,-^ ± )F*(-gi,-g 2 2 ± )(A; + -p + )(fc + -g + ) ejjefc ^g fc 

x (AX a {p)A- A , b {q)Ki k ~ P) A U k ~ «)> ( 21 ) 

where is the Levi-Civitta antisymmetric tensor with i,j = 1,2 and where we defined 
A —A -4- A a 

To obtain the preceding expression we make the assumption that the virtualities in the 
form factors are due solely to transverse momentum such as F(p 2 ,p%) = F(— p\, — p\ ± ). 
This approximation is necessary to recover fc^-factorization in the dilute limit as seen in 
section IIII El 

It is convenient to separate B(k) in four different terms such as 



R f7\ RA 2/L+N2 f d P d( l d 2 p ± d 2 q ± 2 2 \ 77* / 2 2 \ PY 

B z , z >{k) = 64g {k + ) / — — F{-p ± ,-p 2 ± )F (-q ± , - q ^ ± ) €ij€kl — - 

J \ Z71 ) P± ( l± 



dp dq d 2 p ± d 2 q ±r ^ r 2 2 , nf/ 2 2 V\ k 
x (p% a (p±)pA,b(<l±)Ai: a (k -p)A l z , >b (k - q))\ p+=q+=0 (22) 



where z, z' = {p, pA} and where we performed the integration on p + and q + using the delta 
functions in Eq. ffl8l) . These four terms can be evaluated explicitly by substituting the 
solution of gauge fields Eqs. flTTj) . (fT8l) and flT9]) . 

For the first term, it is a straightforward calculation to show that 

BpM = 6 V / ^^> fe V (- f i -pU)f> { -& -4,) 

xe ij e ki 2222 5bc(2n) 2 $ 2 (k± - p± - r ±) 5 ad (2n) 2 5 2 (k ± -q±- s±) 
P±Q± r ± s ± 

x(p* P , d ( r ^PpA s ±))(PA,a(P±)PA,b(q±))- (23) 
The calculation of the second term is similar but requires some more work. By direct 



11 



substitution of the expression of the gauge fields, we have 



dp d 2 p±d 2 q±d 2 r± 



i h I 

p q q 2 

Xe ij e kl 2 J2 Jl ^2 



(2vr) 7 
1 



pi 



2 

(Pt -«0(pa -ie) r\ 



x 



P±T±l2,± Pi ~ ik+e 
(p* A ,a{P±)pA,b{(l±)Kc(P2,± - r±))(p* PiC (r ± )pp jb (q 2 ,±)) 



(2vr) S (p2,±-r ± )5 ac (pX a (p ± )p Aib (q ± ))(p*(r ± )p Pib (q 2> ±)) 



(24) 



The integration on the longitudinal momentum p can be done by looking at the analytical 
structure of the equation. First, we write a part of the integrand as 



1 



(A; — p) 2 — ik + e 
1 



7 

P'2 



(k + — ie)(k — p — ie) 
Pi 



y3 

r 2 . 



2k+ 



k — p~ — 



2k+ 



— ie 



p+=0 



(k + — ie)(k — p — ie) 



-2- 



(25) 



In the complex plane of p~ , the first term in the RHS of Eq. (J25J) has two poles on the 
same side of the real axis and goes like ^pp- when p~ — ► oo. Thus, closing the integration 
contour in the upper-half plane and using the residue theorem, the integration on p~ of this 
term leads to a zero contribution because the contour at infinity has a zero contribution and 
because the contour does not enclose any singularities. For the second term of Eq. ( 1251) . we 
use the principal part (PP) identity = PP^- T i7r5(x). The integration on the principal 
part is zero because the integrand does not depend on p~ and Mm^^ J R R dpPP-^ = 
while the delta function integration is trivial. We finally get 



B p A,p{k) = 32g A J 



d 2 p±d 2 q±d 2 r±d 2 s± 



(27T) 



F(-pi,-plJF*(-ql,-ql 



X€ ^ '2'2 r 2 c 2 (Pp,d( r ±)PpAs±)) 



X 



(p A a(P±)PA, b (q±)}SU^y6 2 (p 2 , ± - r ± )6 ad (2ir) 2 6 2 (q 2 , ± - s ± ) 



+ {pA, a (p±)pA, b (<l±)U: d (P2,± - r ± )}5 bc (27T) 2 6 2 (q 



(26) 



The calculations of B p A )P {k) and B p A, p A(k) are similar. Going through the same steps as for 
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B PtpA (k), we get 



I 



4 , d 2 p ± d 2 q ± d 2 r ± d 2 s ± 2 2 2 2 

32 ^ / TTTTi ^(-P±, -P 2 ,±) F (-?±. -?2,±) 



(27T) 



pi qfc yj gf> 

X6 ^fc' 2,2,2 g 2 (Pp,d(^J-)Pp,c(a±)) 



x 



(PA,a(P±)PA, b (g±))^ c (27r) 2 5 2 (p 2i± - r ± )5 ad (27r) 2 5 2 (g 2 



(27) 



and 



B pA ,pA(k) = 16g 4 



d 2 p±d 2 q±d 2 r ± d 2 s± 



piqk^j gl 

X ¥« 2,2,2 g 2 (Pp,d( r J-)Pp,c(g±)) 



X 



(PAa(P±)PA, b (g±))^ c (27r) 2 5 2 (p 2i± - r ± )5 ad (27r) 2 ( 5 2 (g 2i± - s ± ) 



-(p*A,a(P±)pA,b{ ( l±) U bc{q2,± ~ S ±))5 ad (2n) 2 8 2 {p 2 ,± - r ± ) 

-(p*A,a(P±)pA,b{q±)U* d {p 2 ,± - r ± ))5 bc (27T) 2 5 2 {q 2 ,± - s±) 
+(pla(P±W(?±)^ a *d(P2,± - r ± )U bc (q 2 ,± - s ± )) 



(28) 



The final result for the correlator B(k) = B P:P (k) + B P:pA (k) + B pA:P (k) + B pA , P A(k) can be 
written compactly as 



£(fc) = I6g 4 



4 /" d 2 p±d 2 q±d 2 r±d 2 s± 
(2^)8 



F(-p ± , -p 2 , ± )F*(-q ± , -q 2y± )e tj e kl - 



x (p;,d( r ±)P P ,c(s±)) (PA l0 (p±)pA ) 6(g±)CC(P2 1 ± - rj.)t/ 6c (g 2l ± - s ± )) (29) 
where we defined 



U ab (k ± ) = U ab (k ± ) + (27rY5 2 {k ± )S ab . 



(30) 



From this equation for B(k), we can now evaluate the differential cross section. It is given 
by 

d 2 p± d 2 q± d 2 r ± d 2 s± 



-F(-pi,-pi ± )F*(-qi : -q 2 ) 



Cj^ 



Xe ij e kl' 2 ' 2 2 2 (P*p,c( r ±)PpA s ±)) 

x(PA, a (P±)PA,b(?±)C/* c (p2,± - r ± )U bd (q 2:± - s ± )). 



(31) 
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hnnnr - %ttt + 




+ (-) 




+ (•••) 



pA 




+ (-) + 




+ (-) 

FIG. 1: Diagrams included in the gauge field A A ,Ap and A^ A . The crossed circles tg) represent 
insertions of the strong source PAa( x ±) while the crossed squares M represent insertions of the 
weak source p p , a {x±). By solving the Yang-Mills equation with retarded boundary conditions, it 
resums all the tree diagrams such as the ones depicted in this figure [55]. 



Thus, we can relate the rf production cross section to a correlator of sources and Wilson 
lines. This will be evaluated in the next section in the MV model. Clearly, this expression 
does not have the fc^-factorization structure but it can be recovered in the dilute limit. First, 
we can look at the physical interpretation of the different terms in the cross section. 



D. Diagrammatic Content and Physical Interpretation 

To facilitate the physical interpretation, it is convenient to interpret the gauge field 
expressed in Eqs. (fT7j) . f|T8|) and (fl9l) in terms of Feynman diagrams as shown in Fig. [IJ 
The strong gauge field A A corresponds to a resummation of tree diagrams with any number 
of strong source insertions | |55| . The weak gauge fields A£ and A^ A also resum an infinite 
number of tree diagrams, the difference being that they contain one weak source insertion 
50|. 
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FIG. 2: Diagrams included in rj production at leading order. The thick lines represent insertions of 
the proton (p), nucleus (A) and produced field (pA), and the dashed line is the rj meson. The field 
Ap A contains multi-scattering diagrams shown in Fig. [3j The first figure represents the interaction 
between two gluons producing a rj meson that goes through the nucleus without interacting. The 
second figure corresponds to the situation where the gluons emitted by the proton goes through 
the nucleus before producing the rj . 

P G P 

pLJLJL = ^^JLJL + (...) + 

pA C 

A 

FIG. 3: Multi-scattering diagrams included in A^ A shown in FigO The last diagram is a typical 
diagram associated with the Wilson line. This field contains the multiscattering effects. 

Then, the rj production and the correlator given in Eq. (12TI) can be represented diagram- 
matically in Fig. [2] and [3j These figures show all the diagrams included in the calculation. 
The first term in the figure corresponds to the part where the rj is produced from gluon 
fusion. The second term contains multiscattering effects and eventually, saturation effects. 
Overall, this leads to the following physical picture. A gluon inside the proton interact 
with the classical background field of the nucleus and gets multiscattered. Once it has gone 
through the nucleus, it combines with a gluon and produce the rj . 

E. Recovering ^^-factorization in the cross section 

It is possible to recover fc^-factorization from Eq. (l3Tj) by looking at the dilute limit of 
the nucleus characterized by a weak source such as pA, a <C 1. In that case, we are allowed 
to keep only the first term of the Wilson line expansion U ae (k±) = 2(2n) 2 5 2 (k±)5 ae + O(pa)- 
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In this low-density limit, the cross section becomes 

/o \3o ZJi ^low-density 4 f ,2l f d P-i-d 9-L rp/ 2 2 \ 171* / 2 2 

(27rf2£ fc — = g db ± F(-p ± ,-p 2± )F {-q ± , -q 2 



d 3 k "J (2tt 

i h i 1 

f p C£ < p C£ 

XCij-eW o 2 2 2 2 2 (Pp,a(P2,J_)Pp,6(g2,i_))(PA,a(P±)PA,6(g±))-(32) 

p±y±P2,±y2,± 

In terms of Feynman diagrams, this expression corresponds to neglecting all the multi- 
scattering diagrams shown in Fig. [31 The neglected diagrams are the ones that break 
fcj_-factorization as can be seen from the following argument. 

In this dilute limit, the correlator of the nucleus can also be related to the uPDF 02 like 
in Eq. (jUJ). Using this results, we find that the cross section is given by 

(2nf2E k d(Tl ™ d3 * ensity = j d 2 p ± d 2 q± \F(-p 2 ± ,- q 2 ± )\ 2 

xtijtki — 2~^— <f>i{q±)(f>2{p±)8 2 {k± -p±- q±). (33) 
P±1± 

This is the fcj_-factorized expression of the cross section and is totally equivalent to Eq. 
(IA7p obtained directly from the fc^-factorization formalism (see Appendix EJ. Thus, in 
the low-density limit of pA collisions, we recover a formalism that describes pp collisions in 
;he semihard regime. This is very similar to quark and gluon production in pA collisions 



4a . l50l . |5g, |57|. Note that to obtain this result, it is necessary to assume from the beginning 



that the form factors depend only on transverse momenta. 



IV. COMPUTATION OF CORRELATION FUNCTIONS IN THE MV MODEL 

In this section, we compute the relevant correlators appearing in our expression of the 
cross section using the McLerran-Venugopalan (MV) model. Throughout this calculation, 
we use the notation of 
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59] . We are interested in correlators containing both Wilson lines 
and color charge densities such as the ones included in Eq. (13TI) . In App. [Dj we discuss 
the general case and give more details on the calculation. Note here that to make sense 
of the ordered path, we start with color charge densities that depend on the longitudinal 
coordinate x + . At the end of the calculation, we take p a (x + ,x±) = S(x + )p a (x±) since we 
use the MV model which assumes that the nucleus is moving at the speed of light. In this 
more general case, the 2-point function is simply 

(p a (x + , x ± )p b (y + , y±)} = 5 ab p 2 (x + )5(x + - y + )5 2 (x ± - y ± ) (34) 

where p 2 (x + ) is the average color charge density at point x + . It is related to the average 
color charge density by p 2 = J dx + p 2 (x + ) = A/2nR 2 where R is the radius of the nucleus. 
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In this model, W[p] is still Gaussian, so all even-point functions can be written in terms of 
the 2-point function using Wick theorem and all odd-point functions are zero. 
The Wilson line is defined as 

V 1 



U ab (b + , a + \x±) =V + exp 



dz + / d 2 z±Go(x± — z±)p c (z + , z±)t c 



(35) 



06 



where t c are the SU(N C ) generators in adjoint representation, V + is the path ordering in 
the light-cone coordinate z + and Go is a Green function solution of 

d 1 



-G (x ± ) = 5 2 (x ± ] 



(36) 



With these definitions, it is possible to compute the needed correlators. The general 
strategy is to express the correlators in terms of the following known results for Wilson lines 
in adjoint representation 59(: 



(U ab (b + ,a + \x ± )) = U(b + ,a + \x ± )S, 



Sab exp 



-^-L(x,x)p 2 (b + ,a^ 



(37) 



{U a b(b + , a + \xi±)U cd (b + , a + \x 2 ±)) 



^2 _ ^ V{b ,a \x 1± ,x 2 ±) 



acVbd 



-exp [-N c fc 2 (b + ,a + ) (L(0,0) - L(x 1± ,x 2± ^8) 



where 



L(x,y) = J d 2 z±G (x± - z ± )G (y± - z±) 



(39) 



and where we defined the quantity /j 2 (b + ,a + ) = dz + /j 2 (z + ). In the next subsections, 
we express the correlators appearing in the r( cross section in terms of U and V which are 
defined in Eqs. (|37j) and (|38|). 



A. 1 Wilson line - 1 color charge density correlator: the main building block 

The first correlator does not appear explicitly in the rj' production cross section but it is 
useful to understand the other results. We want to evaluate 

F 1 ' 1 (6+,a+) ee (U ab (b + ,a + \x 1± )p Cl (yf,yi±))- (40) 

where we assume that b + > y\ > a + . For clarity, we define the sources included in Wil- 
son lines as internal sources as opposed to external sources which appear explicitly in the 
correlator (like p Cl (yf,yi±) in Eq. fj40|) ). 
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= 



__^^__^>^_ * o 

FIG. 4: These are the different types of possible contractions using Wick theorem and path ordering. 
The first type is a contraction like (pa j - 1 (zf_ 1 ,z j - 1 ±)pa j+1 (z^ +1 ,z j+1 ±))(p c (y + ,y±)pa j (zf,Zj±)) 
where p c is an external source. The other ones have only internal sources. The second 
one is like (p a . (zj, z j± )p aj+2 (z+ +2 , z j+2 ±)) {p aj+1 (Zj+i, z j+ i±)p aj+3 {z+ +3 , z j+3 ±)) , the third one is 
like (p aj (Zj' , z j± )pa J+3 {Zj +3 , z j+3± ))(p aj+1 (zf +v z j+li _)p aj+2 {zj +2 ,z j+ 2i_)) and the last one is like 
{Pa 3 {zj~, Zj±)p aj+1 {zf +1 , z j+ x±)) {p aj+2 (zf +2 , z j+2 ±)p aj+3 f>+ +3 , z j+3± )). Only the last one has a sup- 
port and thus, a non-zero contribution. 

The first step is to expand the Wilson line. The expression can then be written as 



oo n „ 

F 1 ' 1 (6 + ,a + ) = J2(-9 2 y i Il / ^iGo(xii-^i)(/ fll .../ a „ 

n=0 i=l J 



lab 



X 



'n-1 



b+ 

dzf I dz 2 ... I dz^ 



x(pc 1 (yt^yi±)pa 1 (zt,Z 1± )...p an ( y Z^,Z n± )) (41) 

where we used tl c = —ifabc and where f a bc is the antisymmetric SU(N C ) structure constant. 
Using Wick theorem, the n-point correlation function can be expressed in terms of 2-point 
functions as the sum of all possible contractions. At first, we look only at the contractions 
of Pci(vt,yi±) with p ai (zi ,zx±)...p an (z^,z n ±). Using the MV expression for the two point 
function Eq. (IMl) . we get 



F 1,x (b + , a + ) = ^(yf)G (x 1± -y 1± )J2(-9 2 ) n Yl 

n=0 j=l 
X if ai ••■/a,_i fcxfaj+i ■■■fa n )ab 



n „ 

II / d 2 z i± Go(x 1± - z i± ) 



X 



dzf / dz 2 ... / / dz^ +x ... I dz+ 

a-f J a+ J y\ Ja+ J a+ 

X(p ai (zf, Z 1± )...p a ._ 1 (zf_ 1 , Zj-l±)Paj +1 (Zj+i, Zj+llj-PaniZn, Znl))- ( 42 ) 



As argued in [58|, |59|, only adjacent sources can be contracted due to the ordering in z + . All 



the other nested and overlapping contractions have no support. It can also be shown that the 
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contraction (pa^z^, Zj-i±)p aj+1 (Zj' +1 , z j+ i±)) (where p aj (zf, Zj±) is contracted with the 
external source) does not have support either. These properties are shown diagrammatically 
in Fig. HI They can be used to split the correlator in two parts like 



(Pai( z l, Zl±)~-Paj- 1 (Zj~-l, Zj-l±)Pa j+1 (Zj + l, Z j+1±)--Pa n ( Z n > Z n±)) ~ 
(Ah(A , Ztxj—Paj-i^j-li z j-lA-))(Pa J+1 (z1f +1 , Z j+ i±)...p an (z n , Z n ±)). 

By using these properties and by reorganizing the series, we have 
F 1 ' 1 (6+,a + ) = p 2 (yf)G (x 1± -y 1± )f Cldd , 



(43) 



„2\Z 



oo 

1=0 

x / dz\ 
Jyt -'a; 

oo 
- m=0 



Yl / d 2 z i± G (x 1± - zi 

i=l ^ 



(/ai ---fai )ad 



dz^ ■ ■ ■ 



d 4 (Pax ( Z l, Z l±)~ -Par (4 > Z l±)) 



Vi 



,2\m 



(fbl—fbm) 



d'b 



< I dwf 



dWn 



dw+ (p h (wf,w 1± )...p bm {w+ , w m± ) ) 



44) 



This complicated expression is just a combination of Wilson lines that is given more suc- 
cinctly as 



F l '\b + ,a + ) = /i 2 (^)G (x 1± -i/ 1± )/ Cl ^(f/ ad (6 + ,y 1 + |x 1J „))(^( 2 / 1 + ,a + |x 1± )). (45) 

This can be simplified further by using Eq. f|3T|) and the fact that U is an exponential. We 
can get easily that 



F 1 ' (b + , 



p 2 (yf)G (x 1± - yi±)f Cl abU(b + ,a + \x 1± ). 



(46) 



Specializing to the case of a charge distribution moving at the speed of light we have that 
p a (x + ,x±) = S(x + )p a (x±). Integrating both sides by y + , we get 



(U ab (b + ,a + \x 1± )p Cl (y±)) = p 2 A G (x 1± - y 1± )f Ciab U(b + ,a + \x 1± ). 



(47) 



B. 1 Wilson line - 2 color charges densities correlator 

In this subsection, we compute the correlator 

F 2 '\b + ,a + ) = (U aibl (b + ,a + \x 1± )p Cl (yf,y 1± )p C2 (y^,y 2 ±)). 



(48) 



It is possible to devise diagrammatic rules shown in Fig. that can be used to write 
F 2,1 (b + ,a + ) in terms of known quantities. The rules are discussed in more details and 
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b+ 



<U ab (bta + |xJ) = b 



<U ab (bra + |xJU cd (bta + |yJ> 



<U ab (b:a + |xJP c (ztzJ) = b 



yi^K- + ^jCfTT^.j^h. + (...) 





(...) 



+ ^ J cCX^K^XTb + (...) 



FIG. 5: Diagrammatic rules shown with some of the first topologies they resum. Note that for the 
second correlator (U a b(b + ,a + \x±)U c d(b + ,a + \y±)}, only the ladder-like diagram are allowed since 



58, 



591 ]. Finally, in this notation, the light-cone 



all other topologies have no support and are zero 
coordinates are ordered following the arrow, from left (the smallest) to right (the biggest) so that 
b + > a+. 



ci,y? 



b+ 



c 2 ,yt 



+ 




+ 




FIG. 6: Diagrammatic representation of the 1 Wilson line - 2 color charge densities correlator 
F 2 ' l (b + ,a + ). The first figure corresponds to the contraction of two external sources. The other 
figures represent the cases where the external sources are contracted with internal sources. The 
indices ai, &i, ci, C2, d\, d' 1} da, d' 2 are color indices while a + ,b + ,yf,i/2 are light-cone coordinates. 



generalized to all cases in App. [Dl Using these results, F 2,1 (6 + ,a + ) can be represented 
diagrammatically as in Fig. [61 This includes all possible contractions and topologies that 
need to be resummed. According to the rules of App. [DJ it can be written as 
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F 2A (b + ,a + ) = (p Cl (yfiyi±)pc 2 (y2iy2±))(Ua 1 b 1 ( b+ i a+ \ x i±)) 

+ o(vi - yi)^ 2 {yt)^ 2 (y2)fc 1 d 1 d'Jc 2 d 2 d' 2 Ga{yi± - x 1± )G (y 2± - x 1± ) (49) 

x {U hldl (b + , yf\x 1± )){U d r id2 (yf, y£\xi ± )) {U^ (y£, a + \x 1± )) 
+ 8{yt ~ yi)p 2 {yi)^ 2 {yi)fc 1 d 1 d'Jc 2 d 2 d'.GQ{y ll _ - x 1± )G (y 2± - x 1± ) (50) 
x (U bl d 2 (b + , yt\x 1± )){U d r 2dl (y£, yx\x xx )) (U^ ai (yf, a + \x 1± )). 

This can be simplified further by using the explicit expressions given in Eqs. (JMj) . fl37|) and 
(]46j). We get that 

F 2 '\b + ,a + ) = 5 ClC2 5 aib y(yt)5(yt-yt)5 2 (yi±-y2±)U(b + ,a + \x 1± ) 

+ ii 2 (y£)ii 2 (y£)G (x 1± -yi±)G {x 1± - y 2 ±)U(b + , a + \x 1± ) 

^{fc 1 b 1 df C2 da 1 d{yt -y£) + fc 2 b 1 dfc 1 da 1 0{yt -yt)) (51) 

Considering that the nuclei is moving at the speed of light and integrating on both sides by 
Vi and y£, we get 

(U ab (b + ,a + \x 1± )p cl (y 1± )p C2 (y 2 ±)) = 5 ClC2 S aibl p 2 A S 2 (y 1± - y 2 ±)U(b + , a + \x 1± ) 

+ n\G (xx± ~yi±)G (x 1± - y 2 ±)U(b + ,a + \x 1± ) 

^ (fc±bidfc 2 dai 4" fc 2 hidfcidai ) (52) 

where we defined 

/OO /*0O 
du + / dv + p 2 (u + )p 2 (v + ). (53) 
-oo J u+ 

C. 2 Wilson lines - 2 color charges densities correlator 

In this subsection, we compute the correlator 

F 2 ' 2 (b + ,a + ) = (U aibl (b + ,a + \x 1± )U a2b2 (b + ,a + \x 2± )p Cl (yf,y 1± )p C2 (y£,y 2± )}. (54) 

This calculation is similar to the one of the second correlator F 2,1 using the diagrammatic 
rules. Three typical diagrams of F 2,2 (b + , a + ) out of seven are shown in Fig. [7J The first one 
corresponds to the contraction of the two external sources and the other ones corresponds 
the connected part which consists in all contractions between internal and external sources. 
The first diagram is straightforward to compute. It is given by 

F 2 ' 2 (b + ,a + ) = (U aibl (b + ,a + \x 1± )U a2b2 (b + ,a + \x 2± ))(p Cl (yf,y 1± )p C2 (y£,y 2± )). (55) 

Using Eqs. (|35j) and (|58|) . we get that 

F 2 ' 2 (b + , a + ) = (56) 
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b+ i 
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1 a+ 



+ (...) 



c 2 ,y7 



FIG. 7: These are the first few diagrams included in the calculation of F 2,2 (b + ,a + ). They differ 
essentially by the way the sources are inserted between the blobs. 



The second diagram shown in Fig. [7] is given by 

F 2 ' 2 (6 + ,a + ) = 9(yf -yt)^ 2 (yi)^ 2 (yt)fcid 1 d[fc2d2d' 2 Go(yi±-xi±)G Q (y 2 ±-x l± ) (57) 
x (U bldl (b + , yt\x 1± )U b2f3 (2) (b + , yt\x 2 ±))(U d i id2 (yf, yf\x 1L )U ^) ^ (yf, yf\x 2± )) 
x (Ud^ivt, a+ \ x u.)U$) aa (V2i a+ \%2±))- (58) 

Using Eqs. (1571) and (j55j|. we get 

F 2 ' 2 (6 + ,a + ) = 6(yt -y2)fi 2 (yi)fi 2 (yf)fc 1 d 1 d 1 fc 2 d 2 d 2 G (y 1± - x 1± )G (y 2 ± - x 1± ) (59) 
xV(b + , yf \x 1± , x 2i _)V(yf , yf\x 1± ,x 2± )V(yf, a + \x 1± ,x 2± ) 
= (60) 

which is zero because of the color structure. All the other diagrams included in F 2,2 {b + , a + ) 
can be computed in a similar way. There are five other different ways of inserting the source 
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and they all vanish because of the color structure. The only non-zero term is the first one 
and so we have F 2 ' 2 (b + ,a + ) = F^ ,2 (b + , a + ). Finally, with a nucleus moving at the speed of 
light, we obtain 

(u aibl (b + , a + \x 1± )u a2b2 (b + , a + \x 2 ±)p Cl (yi±)pc 2 (y2±)} = ^ 2 (yi± - y2±) 6ciC2 ^ 1 ^ lb2 

xV(b + ,a + \xi±,x 2 ±). (61) 

This concludes the computation of correlators. We are now in a position to evaluate the 
cross section within the MV model. 

V. NUMERICAL EVALUATION OF THE CROSS SECTION 

A. Proton-proton case 

In this section, we evaluate the cross section numerically in pp collisions. This is done by 
using the expression of the cross section given by Eq. (133|) which can be obtained either from 
the dilute limit of the pA result or from /c^-factorization techniques (see Appendix |Aj). In pp 
collisions, the cross section is related to uPDF that describe the distribution of gluons inside 
each protons. There exist many parametrizations of these distribution functions differing 
mainly in the way the evolution equation is solved. Among the most successful ones are (the 
description of these parametrizations can be found in [qJ 1 ): 

• DIG (Derivative of the Integrated Gluon distribution function) 



CCFM (Catani, Ciafaloni, Fiorani, Marchesini) [8|, [60|, |6l|, [62j, |63|, |64 ] 



KMR (Kimber, Martin, Ryskin) [65 1 



These parametrizations are used to compute the cross section for r( production in pp colli- 
sions and to compare with the result for pA collisions at small saturation scale. 



The final result is obtained by integrating Eq. 



(USD using the VEGAS and the CUHRE 
algorithms implemented in the CUBA package 66j. The number of color is set to N c = 3, 
the center of mass energy to y/s « 200 GeV (RHIC) and the mass of rf to M = 0.957 GeV. 
To make a comparison with pA collisions for a wide range of transverse momentum where the 
MV model is valid (see Fig. (jHJ)), we chose the rapidity y = 1. The results of the numerical 
calculation are shown in Fig. (jSj) where we also present the result for the pA case at small 
saturation scale (Q 2 = 1 GeV) and with the proton described by the parametrization CCFM 
J2003 set 3. 



1 We would like to thanks H. Jung for handing us his FORTRAN routine CAUNIGLU which evaluates nu- 



merically all of these parametrizations. It can be found at http://www.desy.de~.jung/cascade/updf.html 
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FIG. 8: Numerical results of the inclusive differential cross section at rapidity (y = 1) and at 
RHIC energy (y / s=200 GeV). The results for pp collisions are scaled by the number of nucleons A 
to make a comparison with pA. The cross section for pA collisions is evaluated at the saturation 
scale Q s = 1 GeV. 



B. Proton- nucleus case using the MV model 

We now consider rf production for pA collisions for which the cross section is given by 
Eq. ( 13TT) . This formal expression can be simplified by using the results of section ITVl where 
the correlators of Wilson lines are evaluated in the MV model. Note here that in the MV 
model, the nucleus is considered as an infinite source of charge in the transverse plane so 
there are no edge effects. In this kind of description, translation invariance in the transverse 
plane is preserved and therefore the correlators have the following property 

(pA,a( X ±)pA,b(y±) U ce(z±)U c > e >(w ± )) = (p A ,a( X ± ~ ™±)pA,b(y± ~ W±)U ce (z± - W±)U c , e < (0)\G2) 



Then, we can write 
„ o „ do 



g 2 * 2 



(N 2 - 1 

x 



d 2 p±d 2 q±d 2 r± 



X (pA,a {x±)pA,b (V± )Uac{zi_)U bc (w± ) ) 



f p^ Cj^ f^ 1 



e V e kl 2 n 2 r 2 M r ±' X )- ( 63 ) 
P±H± r ± 
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To obtain this expression, we used Eq. Q to convert the average on proton sources to an 
unintegrated distribution function and we Fourier transformed the correlator in the second 
term. We also neglected all terms with only one Wilson line because they are numerically 
very small compared to the other terms. This can be seen as follows. First, using translation 
invariance of the correlator, we see that averages with one Wilson line in the cross section 
are proportional to £7(0) (see Eq. (1621 ). This quantity is small because there is an infrared 
singularity appearing in the argument of the exponential in the following way |59j: 



[7(0) 



cxp 



-^L(0,0)Q s 



exp 



d 2 P± 1 
(2^FpI 



0. 



(64) 



This is not exactly zero because the infrared singularity is regulated by non-perturbative 
effects (confinement) and these induce a cutoff of order Aqcd- The saturation scale in a 

<C 1. Thus, the 



exp 



l QCD J 



large nuclei at small-x satisfies Q 2 ^> Aq CD so that £7(0) 
correlators with one Wilson line can be neglected. 

It is then a straightforward calculation, involving some change of variables, translation 
invariance and Eqs. ( |34l) and ( |6~TI) . to obtain 



(o \3ov da 2 2 2 C f d 2 p±d 2 r ± 2 

(27r) 2Ek ^k = n 9 ^ S± J ^27o^ |F( ^' 



-pi 



p±(p± + a 2 ) 1 



X(p p {r ± ,x) 



(27T)^(p 2 , -r ± )+ / d 2 z ± e 



J(p2,±-r±)-z± 



V(z ± ,0) 



(65) 



where S± = (27r) 2 5 2 (0) is interpreted as the transverse area of the nuclei. We also introduced 
here an infrared regulator given by A in one of the denominator \. This is because gluons 
separated by a distance greater than the size of a nucleon (~ 1 fm) are not correlated 
because of confinement so that A = Aqcd- This is not taken into account explicitly in the 
MV model, leading to infrared divergent quantities [67| which are regulated by adding the 
infrared regulator [56j . This is equivalent to define a gluon distribution function for the MV 
model as 



(j) A (x,p±,Q ) 



2 ^ a s (N 2 -l)A 1 



2tt p 2 ± +A 2 ' 



(66) 



We do not introduce a regulator in the other denominator ^- because as seen in Eq. (1A7I) . 



it is part of the matrix element \Ai\ 2 ~ £ij£ki p ^2 which is well behaved in the infrared. 

In the first term of the cross section, one integration can be done easily with the delta 
function. The remaining integrals will be done numerically. The second term can be sim- 
plified further by doing some integrals analytically as shown in Appendix [Bl The first term 
is very similar to the pp cross section and does not involve any saturation scale dependence 
because it does not include any multiscatterings. The Q s dependence occurs in the second 
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term through the expression of V where we define it as 

Ql = (67) 

1. Kinematical range 

The cross section derived in the previous section is restricted to a certain range of validity 
because we are using the MV model. This model can be used when the radiative corrections 
which goes like a s ln(l/x) are not too large, which is when 0.01 < < 0.1. For rj' 
production, the momentum fraction of gluons in the proton and the nucleus are given by 



x p = \ ±e y , (68) 



XA = (69) 

The kinematic range in terms of the rj' transverse momentum and rapidity where 0.01 < 
xa < 0.1 is depicted in Fig. Ofor RHIC energies. 

It is possible to extend the range of validity to smaller values of xa by using the JIMWLK 



equation 
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45] . This renormalization group equation resums the large radiative 
corrections due to very small x modes. The effect is to change the correlation between 
sources and in general, one looses the Gaussian structure of the weight functional W[p). 
This complicates the computation of Wilson line correlators and is outside the scope of this 
article. 



2. Numerical Results 

We present in this section the numerical results for rf production. We integrate numeri- 
cally Eq. (1B6j) and the first term of Eq. fl65l) using the same values for the parameters as in 
pp collisions. The strong coupling constant appearing in the cross section is evaluated at the 
rf transverse mass scale M\. The unintegrated distribution functions chosen for the protons 
is the CCFM J2003 set 3 because it is very successful in the description of other observ- 
ables like charm and bottom production at the Tevatron The results of the numerical 
integration are shown in Figs. (fTOl) and (fITT) for different values for the rapidity. 

In Fig. (fT2"j) . we show the result for the inverse nuclear modification factor R p a for y = 1 
defined as 

(70) 



We use this because the pp cross section is zero at \k± \ =0. 
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C. Analysis 

For the pp cross section, there is a wide range of variability for the cross section depending 
on the parametrization of the uPDF used. This feature can be used to discriminate between 
the different parametrizations to determine the most accurate one. Thus, by combining this 
analysis with experimental data, the rf production becomes another observable that can 



be utilized to constrain models of uPDF. A similar conclusion was reached in [3l|] where 
the /2-meson production was studied. Note however that the results of the cross section at 
small transverse momentum should be treated carefully. As stated in the appendix |AJ the 
validity of the fc^-factorization approach depends on the presence of a large scale compared 
to the QCD scale given in our case by Q 2 ~ M]_. At very small momentum, we have 
that .i x m 23 for which the semihard inequality An rD C Q 2 < y/s that guarantees 

QCD " ^ 

the accuracy of fc^-factorization is only marginally satisfied. Thus, at small momentum, 
the cross section should be seen as an extrapolation of the fcj_-factorized cross section to a 
regime where /c^-factorization cannot be rigorously proven. 

The pA cross section shows a strong dependence on the saturation scale. As Q s becomes 
smaller, the magnitude of the cross section looks more like the pp cross section. On the 
other hand, the cross section diminishes as Q s becomes larger. This can be understood in 
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FIG. 10: Numerical results of the inclusive differential cross section at midrapidity (y = 0) and 
at RHIC energy (y / s=200 GeV). The first curve is the result of the cross section for pp collisions 
scaled by the number of nucleons A. The other curves are the results for the cross section for pA 
collisions for different values of the saturation scale (Q s =l,4 and 10 GeV). 



the following way. The saturation scale is denned as the momentum at which the probability 
of interaction between different parton cascades is of order one 14], ll5j. The partons having 
a transverse size satisfying 5x± ~ J- > -Jr- will have a very high probability to recombine. 
As the saturation scale is increased, there will be more gluons that will be sensitive to these 
nonlinear effects. The final result is that as more and more gluons have a high probability to 
recombine, the growth of their population will decrease and so is the cross section for a given 
set of parameters. This effect is clearly seen in the cross section computed in our model. 
The fact that the rj' production cross section is sensitive to the saturation scale can be used 
to estimate the numerical value of Q s . The current estimate based on HERA data for deep 
inelastic scattering shows that for RHIC energy, it is given by Q s ~ 1 — 2 GeV 14], ll5| . It 
would be interesting to compare our calculation to experimental data to make an estimation 
of Q s based on rf production and see if it is consistent with the previous estimation. 

The nuclear modification factor shows clearly the saturation effects for 0.4 < \k±\ < 3.4 
because in that range, R p a > 1. At large transverse momentum, it approaches the value 
R p a ~ 0.45. One would expect R p a to be one in that range of momentum because it corre- 
sponds to the regime where there are no saturation effects and where the fcx-factorized cross 
section can be used. This discrepancy can be explained by looking at the approximations 
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FIG. 11: Numerical results of the inclusive differential cross section at rapidity y = 1 and at RHIC 
energy (y / i=200 GeV). The first curve is the result of the cross section for pp collisions scaled by 
the number of nucleons A. The other curves are the results for the cross section for pA collisions 
for different values of the saturation scale (Qs=l,4 and 10 GeV). 



made during the calculation. We neglect two terms in the cross section that are proportional 
to one Wilson line (see Eq. ( 1611) ). In the dilute limit (when \k±\ > Q s ), these terms cannot 
be neglected and have a non zero contribution that would make R p a ~ 1. In that sense, our 
approximation fails at large transverse momentum and our result should be treated carefully 
in that regime. 



VI. CONCLUSION 



In this article, the inclusive cross section for rj production in pp and pA collisions is 
computed. The pp case is analyzed for different reasons. First, it serves as a basis to 
validate our model for rf production against experimental data. This model includes two 
main ingredients. The first one is the effective theory which is at the foundation of our 
analysis. We use in our study a very simplified form of this effective theory where the 



vertex is given by Eq. ([3]). As shown in [19|, |20_|, |2l|, |22j, |23|, |24], |25|], this is a very crude 
approximation of the real vertex, so there are some improvements that can be done in this 
direction in the future. The second ingredient are the uPDF parametrizations. As seen 
in Fig. (jSj), there is still a large variability in the predictions made by different uPDF. 
Therefore, rf production could be used to constrain the models of uPDF once it is compared 
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FIG. 12: Numerical results for the inverse nuclear modification factor at rapidity (y = 1) and 
at RHIC energy ( v / i=200 GeV). We show the results for different values of the saturation scale 
(Q s =l,4 and 10 GeV). 



to experimental data. The same conclusion was obtained in 26J] for the exclusive process 
p+p — > p+p+rj' . Finally, from the theoretical point of view, we studied pp collisions to see if 
the cross section can be obtained as the low density limit of the pA cross section. There are 
now many known examples where this can be seen like gluon production 
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56 



681 ]. quark 



production |48|, |57| and tensor meson production [3jj. We have shown that rf production 
also obeys this property and in that sense, it is a consistency check for the approach used 
for the pA case. 

The pA results for the rf inclusive cross section show some very interesting features. First, 
we show that they are sensitive to the value of the saturation scale. This can also be seen in 
the plot of the nuclear modification factor. This property could be used to make an estimate 
of Q s by comparing with experimental data, which is one of the main goal of our analysis. 
This information is very important for the study of other particle production such as quark 
and gluon production in both pA and AA where saturation effects play an important role. 
Thus, a measurement of rf at RHIC would improve our knowledge of gluon distribution in 
a nucleus. 

Throughout the article, we assumed that the gluon fusion process was the dominant 
mechanism in rf production. There is one other production process that could also be 
important. It is the photon fusion where two off-shell photons emitted by the protons or the 
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nucleus merge to give a rf such as 7* + 7* — > rf + X. This is estimated in |26j| and according 
to this analysis, it should be subdominant in the exclusive production. We assume that this 
holds also in our study, although a careful analysis of this process should be performed. Of 
course, our calculation could be made more accurate by investigating this last issue. 
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APPENDIX A: CROSS SECTION IN pp COLLISIONS 

The calculation of the cross section for rf production in proton-proton collisions at RHIC 
can be performed in the fc^-factorization formalism. This formalism can be used in the 
semihard regime where the factorization scale satisfies the inequality Aq CD C Q 2 C For 
rf, this inequality is only marginally satisfied at small transverse momentum because its mass 
is relatively low. It is recovered at larger transverse momentum, around \k±\ ~ 2 — 3 GeV. 
The starting point is the usual formula of the inclusive cross section in /c^-factorization: 



/„ > o 1 ! da pp ~' ri ' x 9 f 1 dx\ dx2 f d 2 q±d 2 p± , . 9 ^ 9s , . 9 ^ 9 , 



da 9 * 9 *~* r, ' x 

x(2n) 3 2E k — (Al) 



where 



(2n) 3 2E k ^ =^ s \M" ^ \ 2 (2n) A 5\p + q - k) (A2) 

is the high-energy limit of the cross section for off-shell gluons g* to on-shell rf mesons, 
s = X1X2S is the fcj_-factorization flux factor 

11a 

, x\ 2 are momentum fractions of gluons 
and 01,2(^1,2) k±,Q 2 ) are unintegrated gluon distribution functions of proton 1 and 2. The 
unintegrated distribution functions are related to usual parton distribution functions of 
gluons (appearing in collinear factorization ) by 

dk\<\){x, k 2 x , Q 2 ) ss xG{x, Q 2 ) (A3) 

where G(x, Q 2 ) is the usual gluon distribution function in collinear factorization. 

To compute the production cross section of rf mesons, the high energy limit of the lowest 
order matrix element between two off-shell gluons and one on-shell rf has to be calculated. 
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This can be evaluated by using Feynman rules where the vertex, evaluated from the in- 
teraction Lagrangian, is given by Eq. ([2]). Note also that in fc^-factorization, the sum on 
polarization tensors is given by [l|, 0] 

J2^iP)^p) = P -j^ (A4) 

where p\ = {0,p±, 0). The sum on polarizations differs from the usual result because we are 
considering off-shell gluons with a virtuality given by p 2 = —p\. The exact form is due to 
the coupling of gluons to partons through eikonal vertices as well as gauge invariance and 
Ward identities [2J. 

In the center of mass frame, the 4-momenta of partons inside the proton moving in the 
±z direction in Minkowski coordinates can be written, as: 

P= (£ ,0,f) : «=(£o,0,-f) (A5) 
Then, the momenta of gluons in the large energy limit (|p_i_|, \q±\ "C \/s) are simply 

f X ly /s X ly /s\ ( X 2y /s X 2 y/s\ 

v = {—>p±>—) ; * = {—>k-—) ^ 

Using this kinematics, it is possible to compute the high-energy limit of the matrix element 



Ai. This can t 

by (see 



12 



13 



ren be inserted in the expression of the cross section which is finally given 



3l| for more details and similar calculations) 



da pp ~* v ' x 47r 2 
(2n) 3 2E k — = T-^— / d 2 p ± d 2 q ± ct )l {x + ,pl,Q 2 )<p 2 {x.,ql,Q 2 " 



d 3 k (N 2 - 1) 



x5 2 (k ± - P± - q ± )\F(-pl, -gi)| 2 ^f^ (A7) 

p±q± 

where x± = ^=e ±v . This expression can be used to study the phenomenology of rj' produc- 
tion in pp collisions. 



1. Limit of Collinear Factorization 

The procedure to recover collinear factorization cross sections from fc^-factorization is 
well-known |48| and will serve as a consistency check for Eq. (IA7I) . The limit 

IrlIj \q±\ ^0 has to be taken in the matrix elements and the integration on the azimuthal 
angle has to be performed. The last step is to use the relation Eq. (1A3|) to make the last 
integral and relate the unintegrated distributions to the collinear distributions. We obtain 
, pp^r/X 7 r 2 M 2 \F(0 0)\ 2 
{2n)32Ek dk~ = 2s{N 2 V) G i (x' + ,Q 2 )G 2 ( X '_,Q 2 )(27rmk ± ) (A8) 

where x' ± = ^e ±y . This expression corresponds exactly to the well-known result for leading- 
order rf production in pQCD collinear formalism j^"" 
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APPENDIX B: SIMPLIFICATION OF EQ. (f65j) 



In this Appendix, we simplify the second term of the cross section and put it in a form 
that can be evaluated numerically. First, using the definition of V given by Eq. (13*51) . we 
can perform the integrals in the exponent and we get 



V(z±,0) = expl -Q 



v z 2 



f (1-7 + ln 



+ ^ Aq ° D 2^3 ( 1-1:2.3.3 



128 



Aqcdp±| 

7 2 A 2 



(Bl) 



where 7 is Euler constant and where 2-F3 (1, 1; 2, 3, 3 is the generalized hypergeometric 
series. The constant Aqcd appears in this expression as an infrared regulator. 
Then, we consider the second term of Eq. (1651) given by 



W 3 2E k ^ = *WaS ± j d^ ± ^^\F(- P l,-pl ± )\%e kl - 



p^ p T* 



X(p p (r ± , x)e 



p 2 ± (p 2 ± + A 2 )r] 



(B2) 



By letting v ± = k± — p± — r± and by using the identity 

i-2-k 

/ d6e lxcos{6) = 27rJ (a;) 
Jo 

(where Jq(x) is a Bessel function of the first kind), we find that 



(B3) 



(27r)32Efc S = 2 ^VaSx j ^^\F(-pl,-pl ± )\%e kl 



x<j) p [(k-p-v)l,x]T(\v ± \). 



where 



r(KI) 



dzz Jq(\v ±\z)V (z) . 



p l p {k — p — v) 3 {k — p — v) 
p\{p\ + A 2 )(k — p — v)\ 

(B4) 



(B5) 



We let q± = p± + v± and we use polar coordinates to get 



a s {Ml)fi 2 A S_ 



2- 



dpdq / dddQ\F{-p 2 , -p 2 — k + 2pkcos{0))\' 



xpq- 



[qsm(6-e)-ksm(6)f 



(p 2 + A 2 ) [q 2 + k 2 - 2qk cos(6)] 
x<p p [q 2 + k 2 - 2qkcos{Q),x]T{\q 2 +p 2 - 2pqcos{6 - 6)|). 



(B6) 



where a s is the strong coupling constant. This expression is integrated numerically and the 
result is shown in section IV B 21 
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APPENDIX C: COMPUTATION OF EQ. ([31]) IN CO VARIANT GAUGE 



In this Appendix, we compute the cross section by using the solution of the Yang-Mills 
equation in covariant gauge. This serves as a check for Eq. (I3TI) and it shows that our 
result is gauge independent. The solution of the gauge field is given by A^(k) = A^ a {k) + 
AjJJfc) + Ap Aa where A^ a {k) is the field associated with the proton (of 0(p p )), A Aa (k) is 
the field associated with the nuclei (of O(pa) ~ 0(1)) and A^ A [k] is the field produced by 



the collision (of 0(p p p A ) ~ 0(p p )) [50|, |53|. The fields Ap 1 and A^ A are weak and are used 



as small parameters to solve the Yang-Mills equation perturbatively. The field A A is strong 



and satisfies A£, A^ A <C A A . The explicit solution is given by [50|, [53] 

, Pp,a{k 

~k 



AtM = 2ng5(k-)^l (CI) 



PA,a{kl 



A~ Aa (k) = 2ng5(k + )^^- (C2) 

Ai (L\ jg / P P ,b(P±) 

pA > a{ ' ' k 2 + ik+e J (2vr) 2 p\ 

x 



C£(k, P± ) [U ab (k ± - P± ) - (2n) 2 5 2 (k ± -p ± )5 ab \ 
+C»(k, P± ) [V ab {k± -p ± ) - (27r) 2 d 2 (k x -p±)6 ab ]} (C3) 



where V a b(k± — p±) is a Wilson line with a ~ coefficient in the exponential (see [50(]) and 

C+(k, P± ) = 2k + ,Cy(k, P± ) = 2^Cl(k,p ± ) = 2k\ (C4) 

The power counting is very similar to the light-cone gauge. We get that the leading order 
contribution to the field-strength correlator is given by 

BZ'(k) = 64 y d ^F\-pl,-pl A )F{~ q l,- q 2 2 ^,e M 

x ( [q\k - q yA Ab (q)Al b (k - g) - q\k - g)+ A A>b (q) A^k - q)} 

x [p k (k - p) l A A * a (p)A+* a (k - p) -p\k-p) + A2MA l l a {k-p)]). (C5) 

Using the expression of the gauge field, it is a straightforward calculation to show that 
Bpa>(k) — B PtP (k). Just like in light-cone gauge, the other terms require more work. We will 
first look at 

ncov /n RA 3 f d 2 p±d 2 q±dq- 2 2 \ 771/ 2 2 n p k {k-p) 1 

B pA , P {k) = 64^ J — F {- P± , -p 2)± )F(-q ± , -^Je^ — - _ 

x ( [g*(* - g)^^ - g) - g** + ^ 6 (* - g)] 

XpA,fe(g±)PA,a(^)/ ? P,a(^ "Pj) ( C6 ) 
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which is obtained by replacing the gauge fields with their explicit expression. Now, we have 
that 



q\k - q) j A; A>b (k -q)- q l k + ^ pAb {k - q) 



(k - q) 2 + ik+e J (2tt 



d 2 r± Pp, e (r±) 



q l {k — q)ir\ 
k~ — q~ + it 



+ 2q i k + r j 



x 



[U be (k± -q±- r ± ) - (27i) 2 5 2 (k ± - q ± - r ± )5 be ] 



(C7) 



which does not depend on V ab (k±). This is required because this Wilson line does not appear 
in the light-cone gauge calculation. Furthermore, the first term is zero when it is integrated 
on q~ because it has two poles on the side of the real axis (similar to the LC calculation 
in section UlI CI) . The integration on q~ of the second term can be done with the principal 
part identity (see again section [Til Cj) . The result is that B pAp (k) = B P A, p (k). A similar 
calculation can be performed with the two other terms and we find that B cov (k) = B[k). 



APPENDIX D: DIAGRAMMATIC RULES AND n WILSON LINES - m COLOR 

CHARGES DENSITIES CORRELATORS 

In this appendix, we discuss the diagrammatic rules used in Sec. HV] to compute the 
Wilson lines - color charge densities correlators. These techniques can be used in any repre- 
sentation of the SU (N) generators, as long as the Wilson line correlators such as 

FJ(b + ,a + \{a},{b}) = (f/ ai 6 1 (6 + ,a + |a;i±)f/a 262 (& + ,a + |a; 2± )...f/ % ^(6 + ,a + |a; J± )) (Dl) 

are known, which we assume throughout the following discussion. Here, {a}, {b} are the sets 
of color indices defined as {ai, a 2 , aj} and {b%, b 2 , bj} respectively. The most general 
correlator we study here is 

F m ' n {b + ,a + ) = (p Cl (yt,yi±)pc 2 {yZ,y2±)-Pc m {y+,y m ±) 

xU aibl (b + ,a + \x 1± )U a2b2 (b + ,a + \x 2 ±)...U anbn (b + ,a + \x n± )). (D2) 

We want to express this correlator in terms of Wilson line correlators shown in Eq. fIDll) . 
This can be done in a general way. First, we start by defining a number of new quantities 
necessary for the computation. We define a quantity that represents a correlator with a 
number j of color charge densities such as 

G (i,2,...,{fc,fc+i},...j+ij+2) = (pc 1 {yt,yi±)pc 2 {yt,y2±).--pc k _ 1 {yt-nyk-i±)pc k+2 {yt+2iyk+2±)-- 

xpc j+2 {yj+ 2 ,yj+2x))- (D3) 
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In this definition of Gv 12 | fc k+1 y J+lj+2 ), the upper index counts the number of sources 
while the lower index indicates in bracket which sources are missing. We also define H 3,m 
in a similar way by 

H li\...,{k,k+i},---,3+hj+2) = (Pci(yi,yis-)Pc2(y2,y2±)-Pc h _ 1 {yk-vVk-i±)Pc k+2 (yt+2^yk+2x)- 

Xp Cj+2 (^ + +2 ,2/ J+ 2±) 

X^a 1 6 1 (& + ,a + kl±)^a 2 6 2 (& + ,a + |s2±)---^a n 6„(& + ,a + |^n±))conn.(D4) 

where j < m and where the subscript conn, indicates that only the connected part of 
the correlator is considered, which means that all external sources are contracted with an 
internal source (note that external and internal are defined in Sec. IIV Al) . Now, there are 
two possible cases: m can be odd or even. 

1. Even m 

For the case where m is even, we can use Wick theorem to write 

Tpm,n(h+ +\ — r~im TjO,n , \ r rim— 2 fj^< n 

r \p ,a ) - u n -f ^ ^(i,...,i-i,{i},m,..,.?-i,{j}j+i,..,m)^ 

i,j,i<3 



I \ r rim— 4 

~*~ Lr (i,...,i-i,{t},<+i,...j-i, 



i-l,{j},j+l,-,k-l,{k},k+l,...,l-l,{l},l+h-,m) 
i,j,k,l,i<j<k<l 

xH 4 ' n 

* 11 {{l,...,i-l},i,{i+l,...,j-l},j,{j-l,...,k-l},k,{k-l,...,l-l},l,{l+l,...,m}) 

i,j,i<3 

+H m ' n . (D5) 

2. Odd m 

For the case where m is odd, we can use Wick theorem to write 

F m,n (b + , a + ) = ^(l,.. 1 ,i-l,{i},i+l,..., m )-^({l,...,j-l},i,{i+l,...,m}) 

i 

2-^1 (1, {j},i+l,...,i-l,{j},j+l,...,fc-l,{fc},fc+l,...,m) 

i,j,k,i<j<k 

X - n ({l,...,i-l},i,{i+l,...j-l},j,{j-l,...,fc-l},A ; ,{fc-l,...,m}) 
+ ... + ^ G ({l,...,i-l},i,{i+l,...J-l}^0>lv..,m} 

+ J ff m ' n . (D6) 

Then, we have to evaluate explicitly the connected part of the correlator H 3,n . For the 
following, we use the definition 

H m ' n (b + ,a + \{a},{b}) = (p cl {yt,yi±)pc 2 {yt,y2±)-p Cm {y m ,y m ±) 

x U ai bj. (b + , a + 1 x 1± ) U a2b . 2 (b + , a + \ x 2 ± ) . . . U anbn (b + , a + \ x nJ _ ) ) con (iP7) 
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FIG. 13: Splitting the connected correlator H m ' n (b + , a + |{a}, {&}), shown in the first diagram, into 
different slices according to Eq. (|D9|) . Here, aj,bj,Cj,a l j are color indices and a + , b + , cj~ , are 
variables in the + coordinate. 



where we assume for simplicity that b + > y{ > y% > ••• > Vm-i > Vm > a+ - To get the 
general result where all the y^ can have any value between b + and a + , one has to sum over 
all orderings. All other cases where a source is missing such as in the definition Eq. (lD4j) can 
be treated in a similar way. It is now convenient to define a new set of variables which 
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FIG. 14: Inserting a source to compute H 1 ' n (b + , a + \{a}, {b}), which is shown on the left diagrams. 
The (...) contains all the other ways of inserting the source on the other branches. j3 l are color 
indices. 

have the following property y^ > c+ > y^ +1 . Then, we use these new variables to split the 
+ coordinate in slices, which allow us to write the Wilson lines (using their properties) as 



U ajbj {b + ,a + \x j± ) = U v)(b + ,ct\xj±) 



m-2 



U a U) 6 .(c^_i,a + Mp8) 

m + 1 J 



Color charge densities in different slices cannot be contracted using Wick theorem because 
such contractions do not have support. Using this fact and Eq. (ID8j) . we get that 



H m > n (b + ,a + \{a},{b}) = ^(6 + ,c+|{a},{ ai }) 



'm-2 



n^ iin «,c+ +1 iK},K +1 }) 



P =i 



(D9) 



Eqs. (ID7p . (lD9p and the process of slicing the + coordinate are shown diagrammatically in 
Fig. QJ 

The only unknown quantity left is H 1,n (b + , a + \{a}, {b}), to which we now turn. Remem- 
ber that in this quantity, there is only one external source and this source, using Wick 
theorem, can be contracted with an internal source in any of the n Wilson line included in 
H 1,n (b + , a + \{a}, {&}). In Sec. IIV A\ we show that when a source is contracted with a Wilson 
line, it breaks the Wilson line in two, change the color structure and multiply the overall 
expression by ti 2 {y + )Go(x±). For the general case, it is shown in Fig. (|14p . Using this and 
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the definition of H, we get 



n 



H^(b\ci\{a},{b}) = J> 2 (y + )G (y. 



'_L — %j±)fcdd' 



xF n (b + ,y + \{a},{P}) 



F n (y + ,a+\{(3},{b}) 



p j= d> 



. (D10) 



which completes our calculation since the F n, s are assumed to be known from the beginning. 
This general result can be summarized into diagrammatic rules to compute any F m,n . 

• Determine a +-ordering for the external sources. 

• Draw all possible ways the external sources can be contracted with the Wilson lines 
and with each other by respecting the +-ordering. 

• A blob with n external line is given by F^(b + , a + \{a}, {&}). 

• A source insertion at y + with a color index c on line j multiplies the expression by 
[i 2 (y + )Go(y± — Xj±)f c dd' and changes the indices on the left blobs to /3^') = d and on 
the right blob to j3 ij) = d! . 

• External sources that are contracted together gives fi 2 (yi )S ClC2 5(y^ — y2)5 2 {yii. — U2±)- 

• Sum on all +-orderings. 

Examples of these rules are shown in Sec. IIVI 
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